!############################# Change Log ##################################
! 5.0.0
!
!###########################################################################
!  Copyright (C)  1990, 1995, 1999, 2000, 2003 - All Rights Reserved
!  Regional Atmospheric Modeling System - RAMS
!###########################################################################


Module turb_constants

! Nakanishi and Niino (2004) set of constants
real, parameter :: a1=1.18, a2=0.665, b1=24.0, b2=15.0, c1=0.137
real, parameter :: c2=0.70, c3=0.323,  c4=0.0,  c5=0.2
! Constants to define the Turbulent length scale
real, parameter :: q1=3.7, q2=2.7, q3=100., q4=0.2, q5=0.23, q6=5. 
! Constants defined
real :: gama1, gama2, f1, f2, rf1, rf2, ri1, ri2, ri3, rfc
! Some constants to compute E1, E2, etc. from Nakanishi and Niino (2004) equation 20
real :: ce1a,ce1b, ce2,ce3,ce4,cr1
! Part of the constants defined in the Janjic (2001) report, adapted to Nakanishi/Niino
real :: o1, o2, o3, o4, o5, o6, o7, o8, aeh, aem, req,rsl, macheps


!srf- ML -> Constants added for use in nakanishi subroutine
real, parameter ::                    &
        sigwmin     = 1.e-4              &  ! Minimum value of sigma-w
    ,   abslmomin   = 1.e-4              &  ! Minimum abs value of Obukhov length
    ,   ltscalemax  = 1.e5               &  ! Maximum value of Lagrangian timescale
    ,   abswltlmin  = 1.e-4              &  ! Minimum abs value of Theta*
    ,   lturbmin    = 1.e-3                 ! Minimum abs value of turb. lenght


Contains

   subroutine assign_const_nakanishi()

   implicit none

   gama1= 1./3.0-2.*a1/b1
   gama2=(2.*a1*(3.-2.*c2)+b2*(1.-c3))/b1
   f1=b1*(gama1-c1)+2.*a1*(3.-2.*c2)+3.*a2*(1.-c2)*(1.-c5)
   f2=b1*(gama1+gama2)-3.*a1*(1.-c2)
   rf1=b1*(gama1-c1)/f1
   rf2=b1*gama1/f2
   ri1=a2*f2/(2.*a1*f1)
   ri2=rf1/(2.*ri1)
   ri3=(2.*rf2-rf1)/ri1
   rfc=gama1/(gama1+gama2)
!------------------------------------------------------------------------------------------!
! Parameters to compute E1,E2,E3,E4, and R1, from Nakanishi and Niino (2004) equation 20.  !
!------------------------------------------------------------------------------------------!
   ce1a=6.*a1*a1
   ce1b=9.*a1*a2*(1.-c2)
   ce2 =3.*a1*(4.*a1+3.*a2*(1.-c5))*(1.-c2)
   ce3 =6.*a1*a2
   ce4 =12.*a1*a2*(1.-c2)+3.*a2*b2*(1.-c3)
   cr1 =a1*(1.-3*c1)
!------------------------------------------------------------------------------------------!
! New parameters based on Janjic (2001) constraint in L/q for unstable cases               !
! They correspond to the parts of C and D                                                  !
!------------------------------------------------------------------------------------------!
   o1=9.*a1*a2*a2*(12.*a1*(1.-c2)+3.*b2*(1.-c2))
   o2=18.*a1*a1*a2*(b2*(1.-c3)-3.*a2*(1.-c2)*(1.-c5))
   o3=6.*a1*a1
   o4=3.*a2*(7.*a1*(1.-c2)+b2*(1.-c3))
   o5=27.*a1*a2*a2*b2*(1.-c2)*(1.-c3)
   o6=54.*a1*a1*a2*b2*c1*(1.-c3)
   o7=18.*a1*a1*c1
   o8=9.*a1*a2*(1-c2)+3.*a2*b2*(1.-c3)
   aeh=9.*a1*a2*a2*(3.*(1.-c2)*(4.*a1+b2)+a2*b1)
   aem=3.*a1*a2*b1*(3.*a2*(1.-c2)*(1.-c5)+3.*b2*c1*(1.-c3)+6.*a1*c1*(3.-2.*c2)-b2*(1.-c3))+&
       18.*a1*a1*a2*(b2*(1.-c3)-3.*a2*(1.-c2)*(1.-c5))
   req=-aeh/aem
!------------------------------------------------------------------------------------------!
! macheps is the machine precision                                                         !
!------------------------------------------------------------------------------------------!
   macheps=1.
   do while(1. /= 1.+macheps)
     macheps=macheps/2.
   end do
   macheps=sqrt(macheps*2.)
   rsl=(1.+macheps)*((o5+o6*req)/(3.*o1+3.*o2*req))
   
   write (unit=*,fmt='(a)') '--------------- Nakanishi-Niino / Janjic constants ---------------'
   write (unit=*,fmt='(4(a,1x,f12.5,1x))') 'a1=   ',a1,    'a2=   ',a2,   'b1=   ',b1,  'b2=   ',b2
   write (unit=*,fmt='(4(a,1x,f12.5,1x))') 'c1=   ',c1,    'c2=   ',c2,   'c3=   ',c3,  'c5=   ',c5
   write (unit=*,fmt='(3(a,1x,f12.5,1x))') 'q1=   ',q1,    'q2=   ',q2,   'q3=   ',q3
   write (unit=*,fmt='(3(a,1x,f12.5,1x))') 'q4=   ',q4,    'q5=   ',q5,   'q6=   ',q6
   write (unit=*,fmt='(4(a,1x,f12.5,1x))') 'gama1=',gama1, 'gama2=',gama2,'f1=   ',f1,  'f2=   ',f2
   write (unit=*,fmt='(3(a,1x,f12.5,1x))') 'rf1=  ',rf1,   'rf2=  ',rf2,  'rfc=  ',rfc
   write (unit=*,fmt='(3(a,1x,f12.5,1x))') 'ri1=  ',ri1,   'ri2=  ',ri2,  'ri3=  ',ri3
   write (unit=*,fmt='(3(a,1x,f12.5,1x))') 'ce1a= ',ce1a,  'ce1b= ',ce1b, 'ce2=  ',ce2
   write (unit=*,fmt='(3(a,1x,f12.5,1x))') 'ce3=  ',ce3,   'ce4=  ',ce4,  'cr1=  ',cr1
   write (unit=*,fmt='(4(a,1x,f12.5,1x))') 'o1=   ',o1,    'o2=   ',o2,   'o3=   ',o3,  'o4=   ',o4
   write (unit=*,fmt='(4(a,1x,f12.5,1x))') 'o5=   ',o5,    'o6=   ',o6,   'o7=   ',o7,  'o8=   ',o8
   write (unit=*,fmt='(4(a,1x,f12.5,1x))') 'aeh=  ',aeh,   'aem=  ',aem,  'req=  ',req, 'rsl=  ',rsl
   write (unit=*,fmt='(a,1x,es12.5)')   'macheps=',macheps
   write (unit=*,fmt='(a)') '------------------------------------------------------------------'
   
   
   end subroutine assign_const_nakanishi
End Module turb_constants
